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Abstract 

In this article we study a simple, purely topological, cellular model which is allowed to 
evolve through a Glauber-Kawasaki process. We find a non-thermodynamic transition to a 
glassy phase in which the energy (defined as the square of the local cell topological charge) 
fails to reach the equilibrium value below a characteristic temperature which is dependent 
on the cooling rate. We investigate a correlation function which exhibits aging behaviour, 
and follows a master curve in the stationary regime when time is rescaled by a factor of 
the relaxation time tr- This master curve can be fitted by a von Schweidler law in the 
late /3-relaxation regime. The relaxation times can be well-fitted at all temperatures by an 
offset Arrhenius law. A power law can be fitted to an intermediate temperature regime; the 
exponent of the power law and the von Schweidler law roughly agree with the relationship 
predicted by Mode-coupling Theory. By defining a suitable response function, we find that the 
fluctuation-dissipation ratio is held until sometime later than the appearance of the plateaux; 
non-monotonicity of the response is observed after this ratio is broken, a feature which has 
been observed in other models with dynamics involving activated processes. 

1 Introduction 

Over the last few years, much work has taken place on the topic of the glass transition and the 
behaviour of supercooled liquids, both experimentally and through computer simulations. A great 
deal of this activity has involved testing the predictions of Mode-coupling Theory (MCT), which 
has proved thus far to reliably describe many key aspects of the dynamics of glass formers (for a 
review of these, see However, many questions remain, particularly with regard to the extent 
to which one might expect a certain system to satisfy the predictions of MCT. 

The bulk of the numerical studies have taken place on binary Lennard-Jones models, hard 
sphere systems or amorphous silica, all of which involve a number of parameters that can be 
chosen to suit the study e.g. to emulate experimental systems, or to avoid crystallization etc. By 
contrast, the aim of this paper is to investigate a model with as few parameters as possible, in 
order to establish the extent to which the same features and behaviour can be reproduced. The 
model we use is particularly simple in that it is purely topological, involving no length scale at 
all. There is also the added advantage that it shows no tendency to crystallize. Whilst the model 
is the dual of a two-dimensional atomic model, the dynamical driving force is not the same as in 
that case, and there is no true phase transition. 

The lack of spatial co-ordinates within the model invalidates the investigation of many familiar 
functions, such as the van Hove correlation function, the structure factor and the mean squared 
displacement. However, suitable alternatives have been identified which are, in keeping with the 
ethos of this paper, simpler to calculate. 
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Figure 1: A Tl move: cells w and y gain a side, whilst v and x lose a side. 



2 The Model 

The model we use is that introduced by Aste and Sherrington in it is an amorphous two- 
dimensional tiling of cells, with 3 edges incident on a vertex and 2 cells incident on an edge. The 
average number of sides of the cells is constrained to be 6 by the Euler Theorem [|[^. These 
tilings are topologically stable. 

If Hi is the number of sides of cell i, the topological charge qi = 6 — rii measures the deviation 
from the hexagonal configuration i.e. qi is a measure of inhomogeneity. Topological charge is 
conserved under local rearrangements of the cells. 

We define the energy as follows: 

N 

E = - n^f = li2N (1) 

j=i 

where iV is the total number of cells in the system. This quantity is naturally associated with 
the degree of inhomogeneity in the configuration: the ground state (a perfect hexagonal tiling) 
possesses zero energy. 

The system is allowed to evolve through moves as illustrated in Figure |l|, known in the foam 
literature as Tl moves. These are local topological rearrangements exchanging topological charge 
between four cells. Two adjacent cells lose a side each and move apart, and two second-nearest 
neighbours gain a side each and become nearest neighbours. The energy change associated with 
such a move on 4 cells w, w, a;, y (as shown in Figure |l|) with sides n„, riu,, Ux and Uy is: 

AE{n^,ny;n^,nx) = 2(2 + + Uy ~ - n^) (2) 

We use Glauber-Kawasaki dynamics, which allows evolution of the system even at zero temperature 
provided the move decreases the energy or leaves it unchanged. The probability P of performing 
a Tl move is given by: 

P{n^,ny;n^,nx) = — 7^(1 - 5n,,3)(l - '5„,,3)(1 - <5to,y) (3) 

1 + exp [pAh[niu,ny; Uy^rij:)) 

where (3 is the inverse temperature. The first two (5-functions forbid the production of two-sided 
cells, whilst the last (5- function forbids the production of tadpoles (i.e. self- neighbouring cells); 
unless these forbidden formations are present at the start, they will never appear. 

In the following simulations the system consists of = 9900 cells, with periodic boundary 
conditions. We exhibit results of simulations starting from configurations obtained by randomly 
performing IC^N Tl moves on a perfect hexagonal tiling. This is equivalent to running the system 
a,t (3 — Q, and results in an extremely disordered network, with a value for /i2 of approximately 
13. Time is measured in units of N attempted moves. 
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Figure 2: The behaviour of with temperature. The symbols correspond to the following 
waiting times at each temperature decrement of 5T = 0.05: star - lOiV, cross - 50iV, diamond - 
lOOiV, triangle - 500A^, square - 10007V. 



3 Relaxation Dynamics 

We study the temporal dependence of the energy of the system by quenching from a disordered 
network to T = 2.0 and allowing the system to equilibrate at that temperature for IC^N attempted 
moves. Cooling is then carried out by waiting a time t — at each temperature decrement of 
5T = 0.05. The results (averaged over 3 runs) are shown in Figure ^. The solid line superimposed 
on the graph is the equilibrium curve: the computation of this is covered thoroughly in We 
can clearly see characteristic glassy behaviour: the system fails to reach equilibrium even at slow 
cooling rates for temperatures less than T ^ 0.2, and displays the typical strong dependence of 
the energy on the cooling rate, as found at low temperatures in glasses. 

We also study the correlations within the system. The choice of correlation function for this 
particular system is by no means obvious. In previous work on this topological model, a persistence 
function has been studied, which measures the fraction of the total cells which have NOT been 
involved in a Tl move However, this provides only limited information about the system: 
by the very definition of the function, information about a cell is thrown away once it has been 
involved in a move. If the system was in a metastable state in which it performed then undid 
many Tl moves, this would not be revealed. 

Instead we have have chosen to use self-correlation functions: the first of these is 



EL ("^(^)"6)(n.(r + 0-6) 



G(T,r + t) = -2 (4j 



We add in passing that the results for this particular function differ negligibly in a qualitative 
sense from those for an energy-energy self-correlation function such as that used in conjunction 
with the backgammon model §1. 
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(a) C(t) for 13 = 1.0,2.0,2.5,3.0,3.5,4.0 (from left to 
right). These are equilibrium results. 



(b) C(t,t + t) for (from left to right) /? =4.5, 5.0, 
6.0, 7.0. T = 60, OOOAT. These are non-equilibrium 
results. 
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(c) C(t,t + t) for 13 = 4.0. Prom left to right, the 
curves correspond to r = 10^ iV, lO^iV and 10* AT. 



(d) C(t, t + t) for /3 = 
curves correspond to t 



6.0. Prom left to right, the 
= 10^ AT, lO^iV and 10* AT. 



Figure 3: Correlation functions both in and out of equilibrium (each averaged over 3 runs). 



We shall first present the results in equilibrium, where the correlation functions are stationary 
and possess no r dependence. Thus the correlator C(t, r + t) becomes C{t), a function of t 
only. In order to obtain stationarity, the system was allowed to equilibrate for 60, OOOA^ at each 
temperature before measurements were taken, and for 200, OOOiV at /3 = 4. The energy was found 
to have equilibrated in each case and no r dependence was found. For /3 > 4 it became impossible 
to achieve equilibrium within any practical timescale; even if the curves appeared stationary, the 
energy had not equilibrated. 



Figure 3(a) shows the correlation functions for a range of temperatures. At low values of /3 
(high temperatures) the correlator decays directly to zero. As the value of f3 is increased, we see 
a shoulder develops, indicating that two-step relaxation is taking place. If we were able to attain 
stationarity at higher values of beta we would see this shoulder broaden to a plateau. The plateau 



can be seen more clearly in Figure 3(b) , which shows the correlation function for higher values 



of /3, with the waiting time of 60, OOOA^; however, one must keep in mind that this figure shows 



non-equilibrium results. Figures 3(c) and 3(d) show the strong t dependence that exists out of 
equilibrium. 

We can study the correlation function within the stationary regime. The relaxation around the 
shoulder or plateau, and the early stages of the departure from it, are known as the /3-relaxation 
regime. Mode-coupling Theory predicts the behaviour of conventional correlation functions in the 
late part of this regime (i.e. during the plateau and the departure from it) to be governed by a 
von Schweidler law. This states that in stationary conditions the correlator C{t) can be described 
as follows: 

C{t)=f-B{t/U{T))'' (5) 

where / is the plateau height (also known as the non-ergodicity parameter), and both B and b 
are positive constants. All three should be temperature- independent, and in addition b should be 
independent of the choice of correlator. However, this theory is based on atoms interacting via 
a two-body potential and studies conventional correlation functions with explicit position depen- 
dence, such as the density-density correlator (for a recent review see [|l|). Thus it is interesting to 
plot our rather unusual correlation function against rescaled time t/tr to investigate the predictions 
of MCT in our case. 

Following the procedure of Kob 0, we define the relaxation time tr in Equation |^ as the 
time at which C{t) first drops below e~^. Figure || shows the correlation functions for a range of 
beta against rescaled time. All the functions collapse onto a master curve in the region after the 
plateau or shoulder, thus obeying the time-temperature superposition principle of MCT. In the 
late /3-relaxation regime this curve can be fitted by a von Schweidler law with / = 0.86, B = 0.483 
and b = 0.42 as shown on Figure |^. The data fits well to this law, although because of the noise 
the margin of error on the fitting parameters is large. 



Figure 5(a) shows the relaxation time against inverse temperature. This can be fitted extremely 



well by the following offset Arrhenius function: 

U ^ A + Be^/'^ (6) 



where A,B,C are constants. This is plotted on Figure 5(a) , with A = 3.13, i? = 0.0195 and 
C = 2.85. Thus at very low temperatures (high inverse temperature) the system is displaying 
Arrhenius behaviour, which is characteristic of strong glass formers. 

In the high temperature regime, one can fit the Vogel-Fulcher law typical of fragile glass-formers 



i.e. tr = A epl^'^ ~ with To = 0.202 (see Figure 5(b) ), but it does not suit the data as well 
as the offset Arrhenius curve. 

The relaxation time is also subject to predictions by MCT, which states it should behave as: 

OC (T - T^iCTy^ (7) 

close to a critical Mode-coupling temperature T*^*^^, with T*^*^^ defined through this equation 



By performing such a fit on our data, we find that J'^^C'^ = 0.218, as shown on Figure 5(c) , 
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Figure 4: C{t) for f3 = 1.0, 2.0, 2.5, 3.0, 3.5, 4.0 (from right to left) against rescaled time. The solid 
curve is a von Schweidler fit as in equation H with / = 0.86, B — 0.483 and b — 0.42. 



with 7 = 2.86. Neither this law nor the Vogel-Fulcher law fit the data well very close to T*^'^'^, 
although they both fit reasonably well slightly above T*^*-^-^. However, it would seem that whilst 
one can fit both of these functions in certain regions, the best description at all temperatures is 
provided by the offset Arrhenius law. 

Mode-coupling Theory predicts a relationship between the exponent 7 in Equation |^ and the 
exponent b in Equation They are linked by: 



where: 



^ 2a 2b 



r(H-6)2 _ r(l-a) 



r(l + 26) r(l-2a) 



(8) 



(9) 



Using our value of 6 = 0.42, we find a value for 7 of 3.12, which can be compared with the power 
law fit of 7 = 2.86. These two results agree to within 10%, which is reasonable given that the 
noise of the data makes it is very difficult to fit either exponent accurately. 



4 Response Functions 

A clear indication of the presence of aging is the violation of the fluctuation-dissipation theorem 
[^,|^, which we shall briefly review. In equilibrium the response Ra of an observable A to an 
applied conjugate fleld relates to the appropriate time-correlation function Ca =< ^(i)^(O) > as 
follows (in units of kg = 1): 



RA(t) = -T 



-idCAjt) 
dt 



(10) 



Out of equilibrium the relationship must be generalised since the correlation function is now a 
two-time quantity: 



Ra{t, t + T) = -T-^Xa{t, t + T 



dCA{T,t + T) 

dr 



(11) 



with Xa{t, t + r) deflned by the above equation. If A'^(r, t + r) is equal to 1, equation [Tl| reverts 
to the equilibrium case. 
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(a) The logarithm of tr against inverse temperature. The solid curve 
is an offset Arrhenius fit i.e. U = 0.0195 exp(2.85/3) + 3.13 




(b) The logarithm of U against (T - To)~^. 
The straight line is a best-fit Vogel-Pulcher law with 
To = 0.202. 



(c) The logarithm of the inverse relaxation 
time against the logarithm of {T — Tc). The 
straight line is a best-fit power law with Tc = 0.218. 



Figure 5: The behaviour of the relaxation time with temperature. 



In order to investigate autocorrelations, one can apply a field at time r which has a magnitude 
ho, but is randomly positive or negative across the entire system. That field is left switched on, 
and the integrated response G(r, t + t) is measured, where: 

l-t + T 

G{T,t + T)=hoJ RA{t',t + T)dt' (12) 

For many systems it has been observed that for t and r both large, Xa{t, t + t) depends on t, t 
only through the correlator i.e. Xa{t, t + r) = x(C(t, t + t)) . If that is indeed the case, one finds: 

+ f%iC')dC' (13) 



ho 



Under those conditions the slope of a parametric plot of —TG{T,t + T)/ho against the relevant 
correlation function is the interesting quantity, since: 



dC\ ho 



iC) (14) 



One would expect to find a slope of exactly —1 where the fluctuation-dissipation ratio is upheld. 
Where the ratio is broken, the slope gives information about the form of Xa{t, t + r). 

In our particular case, we apply a perturbation /lo X)t=i ^ 6) at time r. is randomly 
assigned to be +1 or —1, and ho is chosen carefully to ensure linear response whilst also obtaining 
a reasonable signal-to-noise ratio. We follow the integrated response G{T,t + T) = X^iLi ^iip-i ~6). 
It is easy to see that the appropriate correlation function for this observable is of the form of 
that given in equation |4[ although one must be careful with normalisation factors. The data is in 
general very noisy and it is necessary to average over at least 8-10 different field configurations to 
be able to reach any conclusions. 

Before presenting results, we draw attention to the unusual nature of our integrated response 
function. At low temperatures almost all the cells within the system are 5, 6 or 7-sided, with the 
vast majority being 6-sided. This means that a large proportion of the system makes no contri- 
bution to the response whatsoever. Furthermore, processes which reduce the energy by turning 
pentagons and heptagons into hexagons in fact reduce the response. This will be discussed in 
more detail in the following paragraphs; in the meantime, it is sufficient to remember that the 
integrated response function used here has some unusual features. 

Figure ^ shows the results for /3 = 4 and 5. These are non-equilibrium results in all cases, 
as the waiting times are not long enough for the system to have equilibrated. The superimposed 
straight lines have slope -1, but are not best fits; they are there to give the eye a comparative 
reference. These plots should be examined from right to left i.e. short times are on the right, 
where C(t, r -|- t) is close to 1, and long times are on the left. We see that for short times the 
fluctuation-dissipation ratio is upheld; the slope is -1 for all values of r, although the intercept 
does of course vary. In fact, the ratio is not broken until C(r, r -I- t) has decreased to ~ 0.6: 
this is considerably lower than the value of G{t,t + t) at which the plateau occurs. Thus the 
fluctuation-dissipation ratio is obeyed for a timcscale larger than that of the onset of the plateau. 
After the ratio has been broken, curves for different r behave differently, tending towards the 
equilibrium straight line of slope — 1 as r increases. The non-monotonicity displayed clearly in 



Figure 5(a) occurs in every case, although at lower temperatures one cannot run long enough to 
see it, as in Figure |6(b)| . Non-monotonicity of the response has been noted in several other models 
and will be discussed briefly in the next section. In the meantime, with regards to this model this 
unusual behaviour can be understood if we now turn to a discussion of the processes dominating 
the evolution of the system. 

At low temperatures, out of equilibrium, moves which increase the energy occur very rarely; 
we can consider the system to evolve by a combination of zero-energy moves and moves which 
reduce the energy. After a short time (~ 10^ A^) the system consists almost entir ely o f 5, 6 and 



7-sided cells. Under these conditions the possible moves are as shown in Figure |7(a)| and |7(b) 
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Figure 6: ""^'^^^'"^ plotted parametrically against C(t, r + t) at different fi and r (averaged over 
10 runs). 



The moves in Figure 7(c) are considered very unlikely due to the low probability of finding four 



non-hexagons clustered together; this is corroborated by the extreme rarity of finding octagons and 



rectangles present. We see from Figure 7(a) that in certain topological arrangements, pentagons 
and heptagons can effectively diffuse freely through the system; in particular, pentagons and 
heptagons can 'pair up' to move through the system together with no energy cost. To reduce 



the energy, one must annihilate two or more 'defects' (i.e. non- hexagons) as in Figure 7(b); this 
requires the absorption of a 5-7 pair by an isolated pentagon or hexagon, or the annihilation of 
two 5-7 pairs. This brings us to a conceptual picture containing both fast and slow dynamics: 
the fast dynamics is due to the rapid diffusion of 5-7 pairs moving freely through the system, 
whereas the slow dynamics is due to annihilation processes. The latter can be broken down into 
two different types: absorption of a 5-7 pair by either a heptagon or a pentagon (which can also 
rearrange the network such that isolated pentagons and heptagons become 5-7 pairs), and the 
complete annihilation of two 5-7 pairs. Of course, in equilibrium an equal number of 5-7 pairs are 
created as annihilated/absorbed, but in the relaxation to equilibrium, annihilation has the upper 
hand as the system starts from a configuration in which there are more 5-7 pairs present. 

The effect of the field is to try and 'pin' defects onto the appropriate cell i.e pentagons try 
to settle on sites with = +1 and heptagons on e.; ~ —1. Thus with the field switched on, the 
diffusion processes move the defects around until they are on the appropriate sites, and in this 
way increase the value of the response — TG(r, r -I- t)/ho. However, annihilation processes remove 
defects from the system entirely, thus reducing the response and also reducing C(t, r + t). We 
have competition between the annihilation processes and the diffusive processes; when annihilation 
dominates, the slope on the parametric plot becomes positive, as seen for small values of C(t, r+t) 



in Figure 6(a) . If one looks instead at a quantity such as ~TG{t,t + t)/hofi2 (in some sense a 
'response per defective cell') one finds that the slope remains negative: this tells us that the 
diffusion processes are succeeding in placing more of the defects that remain on the appropriate 
sites, but the number of those defects present in the system is declining. The value of the response 
reached when the system finally achieves equilibrium depends on the waiting time r, and increases 
as r increases; this is because both the response and the correlation function have been normalised 
by a factor that is equivalent to the energy at time r (see equation ^) . 

This picture of diffusion and annihilation/creation processes is helpful to an understanding of 
the plateaux in the correlation functions at low temperatures. The initial descent to the plateau is 
due to diffusing 5-7 pairs, which quickly move the network away from the starting configuration. 
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(a) Zero energy Tl moves. (b) Tl moves which reduce the en- (c) UnUkcly zero energy Tl 

ergy. moves. 

Figure 7: The possible Tl moves at low temperatures. 

However, the lone defects are trapped, and can only be freed by a move which costs energy or 
by absorption of a 5-7 pair; both of these processes occur on a timescale which is temperature- 
dependent, and thus the length of the plateau itself is dependent on the temperature. This 
picture of trapped or caged defects is conceptually similar to that typically used when dealing 
with Lennard- Jones binary models Q pO[ |. 

5 Concluding remarks 

In summary, we have studied a simple two-dimensional topological glassy model, testing some 
of the predictions of Mode-coupling theory for this system. We find that this theory can indeed 
be used to describe the behaviour of the system in so far as we have investigated, although this 
description fits some features more closely than others. The correlation function follow a von 
Schweidler law as predicted in the late /3-relaxation regime, and the exponent of this roughly 
agrees with the exponent of the power law fitted to the relaxation time data. However, the 
relaxation times fit better to an offset Arrhenius law in all temperature regimes. 

Investigation of the fluctuation-dissipation ratio reveals an interesting feature, namely non- 
monotonicity of the response. This is a feature which has also been observed in other simple 
models; for example, vibrated granular media pT| , constrained Ising chains and the Backgammon 
model in one dimension p2[ , and also recently in a two-dimensional short-range spin model with 
uniform ferromagnetic 3-body interactions ]l3[ |. In each case, the dynamics at low temperatures 
can be considered to involve activation over energy barriers; as mentioned in ]l3[ |, this raises the 
question as to whether this is a generic feature of models with activated dynamics. 
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